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Abstract 

We study the method for generating the initial data of black hole systems in Gauss-Bonnet (GB) 
gravity. The initial data are assumed to be momentarily static and conformally flat. Although the 
equation for the conformal factor is highly nonlinear, it is successfully solved by numerical relaxation 
for one-black-hole and two-black-hole systems. The common apparent horizon is studied in the 
two-black-hole initial data, and the result suggests that the Penrose inequalities are satisfied in 
this system. This is the first step for simulating black hole collisions in higher-curvature theories. 
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I. INTRODUCTION 



Higher-dimensional gravity has been attracting a lot of attentions motivated by the TeV- 
gravity scenarios jl-3|. If the three-dimensional space is a 3-brane in large or warped extra 
dimensions, the Planck energy could be of O(TeV), and the trans-Planckian collision will 
happen at accelerators such as the Large Hadron Collider (LHC) 4-f|. If this is the case, 
the intermediate state of a collision with sufficiently high energy and small impact parameter 
is expected to be a mini black hole, and this motivated a lot of works (see jz] for a review). 
At this time, no evidence for black hole signals has been found, and the restrictions for the 
Planck energy and the minimum black hole mass Mg™ n ' ) (such that the produced object can 



be regarded black hole if M > M^ n) ) have been derived [8||. Another motivation for 
studying higher-dimensional gravity is the AdS/CFT correspondence, which conjectures the 
duality between gravity in anti-de Sitter (AdS) spacetime and the conformal field theory 
(CFT) on the boundary of that spacetime. If this conjecture is correct, the phenomena in 
CFT for which direct calculation is difficult due to strong coupling effect can be predicted 
by calculating the dual gravitational system. 

One of the important approaches for higher-dimensional gravity is to explore the nonlinear 
dynamics of higher-dimensional spacetimes by numerical relativity. Several formulations 



ft 



Hi, 



and codes of higher- dimensional numerical relativity have been developed so far 
and interesting simulations have been performed: the time evolution of Gregory-Laflamme 

141 , dynamics of complex scalar 



instability 



12J, slow- velocity collision of black holes 
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field minimally coupled to gravity in Kaluza-Klein spacetimes [11] . and bar- mode instability 



of rapidly rotating Myers- Perry black holes 15|, |l6). In these works, time evolutions of 



3S 111, 



spacetimes are simulated in the framework of general relativity (GR), which is the simplest 
theory of higher- dimensional gravity. 

One of the interesting extensions of higher-dimensional numerical relativity is to include 
the higher-curvature terms. In four dimensions, the Lagrangian density that leads to the 
second-order equation of the metric is just the Ricci scalar. However, in higher- dimensions, 
the Lagrangian density including higher-curvature terms also leads to the second-order equa- 
tion if the combination of higher- curvature terms are chosen appropriately. These theories 
are Gauss-Bonnet (GB) gravity or, more generally, Lovelock gravity [isj]. Among these 
higher-curvature theories, we take attention to GB gravity in this paper. The GB term gives 
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a quadratic correction with respect to the Riemann tensor to GR, and this GB correction 



arises in the low-energy limit of the heterotic string theory [19N21|. Therefore, the GB term 
may have important effects on the mini black hole phenomena at accelerators. Also, the 
AdS/CFT correspondence in the context of GB gravity has been considered (e.g., 22l~l24|). 



The black holes in GB gravity have interesting properties. Although some part of GB 
gravity resembles GR, the other part does not. For example, the solution of a static spheri- 
cally symmetric black hole found in Ref. 25( has two branches: the non-GR branch that is 
asymptotically AdS and the GR branch that is asymptotically flat (see H] for a detailec 
study on causal structures of these solution). The GB version of Birkoff's theorem 28|, [29 1 
(see also [30|] for Lovelock gravity) states that the spherically-symmetric vacuum spacetime 
(at least locally) corresponds to one of the two branches. The black hole of GR branch is 
shown to be unstable for D = 5 against scalar mode and for D = 6 against tensor mode if the 
coupling constant «gb of the higher-curvature term is sufficiently larger than appropriate 
power of GM, where G is the gravitational constant and M is the Arnowitt-Deser-Misner 
(ADM) mass of the black hole 31 34 ] (see also j^] for a study that proves the generic 
appearance of instability of a spherically-symmetric black hole in Lovelock gravity). 

The GB version of numerical relativity (say, numerical GB gravity), if it is developed, 
would play important roles to explore nonlinear phenomena in GB gravity. For example, 
the temporal evolution of the instability of a small spherically symmetric black hole could 
be followed by simulations to clarify the final end state. Also, by simulating the time 
evolution from initial data of a nonstationary rotating black hole, it would be possible to 
obtain an indication for the stationary rotating black holes whose analytic solutions have 



not been found to date (but see 36] for a numerical construction of the black hole with two 
equal rotational parameters). Of course, the simulation of high- velocity collision of black 
holes is an interesting issue to clarify the effect of higher- curvature terms on mini black 
hole phenomena at accelerators (see j^] for a partial result for a head-on collision of the 
Aichelburg-Sexl particles in GB gravity). 

There are few studies for numerical GB gravity. The N + 1 formalism, which is the 
extension of the ADM formalism to GB gravity, has been done by Torii and Shinkai 38]. 
Similarly to the ADM formalism, Einstein-Gauss-Bonnet equation is decomposed into the 
Hamiltonian and momentum constraints and the evolution equations. On the other hand, 
the numerically stable formalism of GB gravity which is analogous to, e.g., Baumgarte- 
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Shapiro-Shibata-Nakamura (BSSN) formalism [39|, |40( has not been developed yet, and also, 
no simulation of numerical GB gravity has been done. 

When one simulates a spacetime in GB gravity, similarly to the case of GR, the first step 
is to prepare an initial spacelike hypersurface (i.e., initial data) by solving the constraint 
equations. In this paper, we focus attention to the method for preparing initial data of black 
hole systems in GB gravity. As a first step, the initial data are assumed to be momentarily 
static (i.e., time symmetric). In GR, the Brill-Lindquist initial data [41] are well known and 
widely used as the initial data of momentarily static multi black holes (see also 42|, |43|] for 
higher-dimensional studies). We discuss the extension of the Brill-Lindquist initial data to 
GB gravity, and successfully generate the initial data for one black hole and two equal-mass 
black holes. 



Using the two-black-hole initial data, we discuss whether the Penrose inequality [44J holds 
in this system. The Penrose inequality states that the area of an apparent horizon (AH) 
is not greater than that of the Schwarzschild-Tangherlini black hole with the same ADM 
mass. In the case of the momentarily static initial data, the Penrose inequality was proved 
for 4 < D < 7 using the conformal flow method 45) . It is of interest whether such a bound 
on the AH area holds or not in GB gravity. In addition to the original Penrose inequality, 
we also discuss whether the AH area is bounded from above by that of the spherically- 
symmetric black hole in GB gravity. The results suggest that both of the inequalities are 
held in this system. 

This paper is organized as follows. In the next section, we review GB gravity and N + 1 
formalism in this theory focusing attention to the part that is related to our study. In 
Sec. Ill, general framework for constructing the GB version of the Brill-Lindquist initial 
data is explained. The equation for the conformal factor \1/ is a Laplace equation with 
formal source term which is highly nonlinear in and we prove the regularity of the source 
term which is necessary for the existence of the solution. In Sec. IV, the one-black-hole 
initial data are constructed numerically, and we show the agreement of the data with the 
time-symmetric slice of a spherically symmetric black hole spacetime. In Sec. V, the two- 
black-hole initial data are constructed. Using those data, we analyze the condition for the 
AH formation, and the Penrose inequalities in this system are discussed. Sec. VI is devoted 
to summary and discussion. Throughout this paper, we use the unit c = 1, while the 
gravitational constant G is explicitly written. 
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II. CONSTRAINT EQUATIONS IN GAUSS-BONNET GRAVITY 



In this section, we review GB gravity and the N + 1 formalism of Ref. 38 1 focusing 
attention to the part that is related to the setup of our study. 



A. Gauss-Bonnet action and equations 



In this paper, D denotes the dimensionality of the spacetime Ai with the metric g^. 
We also introduce N = D — 1, which is the dimensionality of a spacelike hypersurface S 
in the spacetime Ai. Throughout this paper, we assume the spacetime to be vacuum. The 
Einstein- Gauss-Bonnet action 



171 | in .D-dimensional spacetime is 



S 



1 
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(TZ + a GB C GB ) y/^gd D x 



M 
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n 2 - wijRr + n uuoa n» up ° 



(i) 



(2) 



where 1Z, TZ^ U , and TZ^ up(T are the Ricci scalar, the Ricci tensor, and the Riemann tensor, 
respectively, and «gb is the coupling constant with dimensionality of squared length. The 
action ([T]) gives the gravitational equation as 



+ ttGB"H 



0. 



where 



and 



K 



(3) 



(4) 



(5) 



In the case D = 4, the term £gb in the action ([I]) gives a topological invariant as proved 
in the generalized Gauss-Bonnet theorem |46|. and hence = 0. The GB term becomes 
nontrivial only for higher dimensions, D > 5. 



B. The N + 1 formalism 



Here, we review the initial value equations in GB gravity based on the N + 1 formalism 
developed in Ref. [38J. We consider the spacelike hypersurface E with the induced metric 



7|ti/ an d the extrinsic curvature K^ u , and let n' 1 be the future-directed timelike unit normal 
to the hypersurface. The projections of the gravitational equation (Q^ u + aGB^^nV ', 
(Qfiu + c ( GB'HfMu) n ^l l 'p, (Gfj,u + c e GB'Hfj,u)l^pl u (T give the Hamiltonian constraint, the momentum 
constraint, and the evolution equations, respectively. The initial data should be prepared 
so that they satisfy the two constraints. The Hamiltonian constraint is written as 



M + « GB (M 2 - AM ab M ab + M abcd M abcd ) = 0, (6) 



where 



Mijki = Rijki + {K ik Kji - KuK jk ), (7) 
and Mjj —_ and M = 7 ab M a b. The explicit form of the other equations can be found 



in Ref . 



38|. 



C. Conformal approach 

Hereafter, we focus our attention to the momentarily static case, or in other words, the 
time-symmetric case: fQj = 0. In this case, the momentum constraint is trivially satisfied, 
and MijM = Rijki- We further assume the initial space £ to be conformally flat: 



^4/(7V-2) 



7ij) 



(8) 



where 7^ is the flat-space metric. Note that the authors of Ref. 38j treated 7^' as an 
arbitrary metric, and here we chose the very special case. Also, in Ref. [38J, the conformal 
factor was set to be \l/ 2m , where m is an arbitrary number. Here m = 2/(N — 2) is chosen, 
because in this case, we can treat the problem as a natural generalization of the GR studies. 
In this setup, the equation for the conformal factor is written as 



D a D a m = a GB S, 



(9) 



where D a is the covariant derivative with respect to the flat-space metric. Here, S is written 



as 



S 



N-3 



N+2 



%-*r=s J 4(jv - 2) (D a D a ^y - (D a D b m)(D a D b m) 



{N -1){N -2) 
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III. METHOD OF INITIAL DATA CONSTRUCTION 



In this section, we give the formulation for generating the initial data with iV B H black 
holes in GB gravity. 



A. Decomposition of conformal factor 

Since we would like to obtain the initial data of GB gravity as continuous extension of 
those of the GR case, we decompose \1/ as 



the deviation from the GR case in the presence of acB ^ 0. The equation for g becomes 



* = + a GB g. 



(11) 



Here, \l/o is the solution in the GR case (i.e., D a D a ^ = 0), and the term cxgb9 represents 



D a D a g = S 



(12) 



with 



N-3 



4 



s = 



(N - l)(N-2) 




(n) 



(13) 



where 



s (o) 



4(7V - 2)(D a D b ^o)(D a D b ^ ) + 8N*-\D a * )(Db*o)(D a D b *o) 

N-2 



y- 2 (m ) 4 , (14) 



8(N - 2){b a b b m Q ){b a b b g) - 8^\bv ) 2 (b a b a g ) 

+ \(b a v )(b b v )(b a b b g) + 2(b a v )(b b g)(b a b b y ) 




^-y- 2 (b* ) 2 (b a * )(b a g), (15) 
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s (3) = -8^-\Dg)\D a b a g) + 8N^-\D a g)(D b g)(D a D b g) 

] *-\D a %)(D*g)(Dg) 2 , (17) 



16iV(Ar-l)_- w ~ ... , a 



iV-2 

and 

, w = _ 4Ar(*-i) -_ a ,^ 



a W = ^3^* ™ (18) 

We call the right hand side S of Eq. (1121) the "source term" hereafter. 

B. Puncture initial data 

Here, we specify the GR solution ty of a system with TVbh black holes. For this purpose, 
we introduce the Cartesian coordinates (x a ) to the flat space and specify Nbh points at which 
\l/o diverges (i.e., punctures). The location of the n-th puncture is denoted by x a = x® n ->, 



X (n) 



and we define x a ^ = x a — % and 
D a D a ^Q = is adopted as 

Nbh 



Then, the solution to the Laplace equation 

H 

# = 1 + ^W), ( 19 ) 



n=l 

with 

= (N- m Z Kr (20) 

where £In-i denotes the area of a (N — 1) -dimensional unit sphere and is a mass 

parameter for n-th black hole. The Arnowitt-Deser-Misner (ADM) mass Mq is given as 



M = J2n=i M o n) - This solution was studied in Ref. j4jj in the D = 4 case and is often called 
the Brill-Lindquist initial data (see 



42|,|43( for the studies in higher- dimensional cases). This 
space possesses the N Bii Einstein- Rosen bridges (say, throats) and each puncture corresponds 
to the asymptotically flat region beyond each throat. 

C. Finiteness of the source term 

The conformal factor ty of the Brill-Lindquist initial data in the GR case diverges at 
each puncture as Rr„^ ^ '■ Then, a naive estimate gives a severely divergent behavior of the 
source term S ~ at each puncture. If this is the case, g also has to diverge at each 

puncture, and then, the behavior of S would be further modified, implying further stronger 
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divergence of g. This procedure would continue eternally, and hence the solution g would 
not exist. 

However, this naive estimate is not correct, since the cancellation of divergent terms 
occurs. As a result, S becomes zero at the punctures, and hence the source term S is well 
behaved. Let us check this in the following. 

Assuming a regular behavior of g at each puncture, the functions s^ 2 -*, s^ 3 \ and obvi- 
ously become zero at each puncture after multiplying the factor i|/~( 7V + 2 )/( iV_2 ) in Eq. (TT3]) . 
Therefore, we focus our attention only to and s^. We substitute Eq. ffl9|) into Eqs. ( I14p 
and (1151) . rewrite it with R{n)4>'L) + (N — l)V'(n) = 0> anc ^ collect only terms for which di- 
vergence at the punctures is suspected. For example, ^^(D^q) 4 from the third term of 
Eq. ffT4"|) is calculated as 

k,Lm,n 



fc k,l(k^l) 



(21) 



where s = xf^/R^) is the unit vector and (ri(fc) • n^) is the inner product. Here, 0(ip' 2 ) 
are the terms that become zero at the punctures after multiplying the factor \[/ _ ( iV+2 )/( 7V " 2 ) in 
Eq. fTl3|) . and thus, we do not consider these terms because we are interested in cancellation 
of divergent terms. In this manner, s^ -* and s*- 1 -* are calculated as 

AN(N - 1) ^ ( $ fe - 



N -2 



(J0_ 



[%)v; fe) + (iv-2)*]- 



k,l(k^l) 



R{k)R{i) 

+ 0(R^), (22) 



-8£[ity)Vfo + (tf-2)¥] 



(k) 



x 



D a D a g + Nn a {k) n b {k) D b D a g + 2A ( A " ' ' ' ' ^ 



iV-2 



Now, we use the relation 



i2 (n) V; n) + (iV - 2)V ( n) = o. 



+o(i^5). 



(23) 



(24) 
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and find = 0(R7^) and = O(R^). Therefore, S behaves as ~ ^nj 2 a ^ eacn 
puncture, and no divergence of the source term occurs. For this reason, we can expect the 
existence of the solution of g that is regular at each puncture. 

In Sees. IV and V, we explicitly construct the numerical solution of g for one black hole 
and two black holes, respectively, assuming the regularity of g at each puncture. There, 
it turns out that this numerical method works well. It is worth pointing out the relation 
between our method and the method for generating non-time-symmetric initial data of 
puncture-type boosted black holes in GR developed by Brandt and Briigmann {47]. In 
their formalism, the space is assumed to be conformally flat, and the analytic solution 
of the extrinsic curvature found by Bowen and York |48|] is used. The conformal factor is 
decomposed as \I/ = ^>q+i/j, where \l/o is the conformal factor for the Brill-Lindquist solution. 
The equation for if) becomes the Laplace equation with a formal source term depending on 
\P. Here, the source term is shown to become zero at each puncture, and hence, ip can be 
solved numerically assuming the regularity at each puncture. Therefore, our method is very 
analogous to the Brandt-Briigmann formalism. 

D. ADM mass 

Here, we discuss how to calculate the total gravitational energy. The ADM mass M is 
given by 

M = -4$r J kL t >°* dS "' < 25) 

where S is the surface at infinity. Suppose \I/ = \I/o is the conformal factor with the ADM 
mass Mo in the GR case. Then, using the Gauss law and Eq. (I12p . we find that the ADM 
mass in the GB case is 

M = M - [ N ~^ a ^ I Sdxt-.JxN, (26) 



4vr(A - 2)G J± 

where £ is the whole flat space. It is important to point out that the ADM mass M in the 
case acB 7^ is different from Mo, although Mo is the ADM mass in the GR case ogb = 0. 
Therefore the parameter Mo is (say) an artificial mass for o>gb 7^ 0, and the true ADM 
mass M is determined after \1/ is solved. The similar phenomena can be found also in the 
Brandt-Briigmann formalism. 
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IV. ONE-BLACK-HOLE INITIAL DATA 



Let us begin our numerical analysis with the one-black-hole initial data. Here we assume 
the initial data to be spherically symmetric and introduce the radial coordinate R in which 
the metric of the flat space becomes 



ds 2 = dR 2 + R 2 dn 2 N ^. 



(27) 



Here, dfl 2 N _ l is the line element of the (N — 1) -dimensional unit sphere. We impose g = g{R) 
and set the function \l>o in Eq. (ITT]) to be 

'Rs(M ) 



1 + 



where we defined Rs(M ) as 



Rs(M ) :-- 



R 



AtiGMq 



N-2 



(28) 



l/(N-2) 



(29) 



Y(N-i)n N . 

Note that Rs{M ) is different from the Schwarzschild radius rs{M ), and they are related 
as r s {M,)=A^ N - 2 )R s {M Q ). 

In the case of «gb — 0, the conformal factor $ = \l/ gives the spacelike hypersurface 
which is known as the Einstein- Rosen bridge (i.e., the time-symmetric Cauchy slice in the 
Schwarzschild- Tangherlini spacetime). The coordinate R is called the isotropic coordinate, 
and the minimal surface (or the apparent / event horizon) is located at R = Rs(Mq). Because 



of the GB version of Birkoff 's theorem 



28 



291 ]. also in the case of ogb > 0, the initial data 



is expected to give the time-symmetric Cauchy slice of the spherically-symmetric black hole 



spacetime 



251 ] whose metric is 



dr 2 



ds 2 = -f(r)dt 2 + + rW 



f(r) 



JV-D 



with 



/(0 



1 + 



2a 



1 + 



GB 



(30) 



(31) 



where «gb = {N — 2)(N — 3)«gb and M = [r s (M)] N ~ 2 . If we can find the coordinate 
transformation from the Schwarzschild-like coordinate r to the isotropic coordinate R, the 
conformal factor ^{R) is obtained. However, since f(r) has a complicated form in the cases 
cxqb 7^ 0, it seems impossible to find the analytic formula for the coordinate transformation 
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in contrast to the GR case. Therefore, this problem has to be solved numerically. Besides, 
we can obtain a lot of lessons for the numerical method from this problem as we will see 
later. 



A. Perturbative analysis 

In the numerical calculation, Rs(M ) is adopted as the unit of the length, and also in 
the following, we use the unit Rs(M ) = 1 unless explicitly specified. Let us consider the 
case where the coupling constant «gb is small, and thus, the correction from the GR case 
can be treated as a perturbation. In this case, $ in Eq. ffT2~]) should be replaced by ^o, and 
the equation is reduced to 

M — 1 Q 4 

9,rr + ~^-9,R = -4iV(iV - 3) R- 2 ^l R % ~ N ~ 2 . (32) 

Here, the cancellation of divergent terms has occurred in the right-hand side because of the 
relation ^ 0tR + (N - 2)^ /R = (N - 2)/R, and as a result, the source term is 0(R N ~ 2 ) in 
the neighborhood of the puncture. In the case of D = 5 (i.e., N — 4), this equation can be 
solved analytically: 



9 = - ,K v.:.;s (33) 



4 (1 + 3R 2 + R 4 ) 
3(1 + R 2 f 

Thus we explicitly confirm the existence of the solution that is regular at the origin R = 0. 
This result leads to the expectation that a regular solution of g exists in more general cases 
including higher-order terms in acB- 

The function g behaves as ~ (4/3)_R -2 at the distant region, R ^> 1. This indicates that 
the ADM mass is shifted as 

M = M + (34) 

Gr 

Therefore, we confirm the statement of Sec. IIII Dl explicitly: If we impose the regularity of 
g at the origin, the function g also contributes to the mass, and hence, the ADM mass is 
determined after g is solved. 

B. Numerical approach 

Now, we study the numerical generation of g(R) for the general cases of ogb > 0. We 
write down the Laplace operator D a D a and the source term S in terms of the radial coor- 
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dinate R. Here, we adopt the formulas of and for which the divergent terms are 
canceled out by the relation ^ 0jR + (N -2)V/R= (N — 2)(1 + a GB g)/R. Then, the right 
hand side becomes 0(R N ~ 2 ) similarly to the perturbative case and the source term S is 
well behaved at the origin R = 0. This confirms the general proof for the finiteness of S of 
Sec. IfflCl 

The numerical calculation is done in a finite region, and hence, there is an outer boundary 
of the computation domain, R = -R max . Here, we have to specify the outer boundary 
condition. As discussed in the previous subsection, the behavior of g at R ^> 1 is expected 
to be g ~ C / R N ~ 2 . Since the value of C is not known before solving g, we eliminate C using 
the combination of g t R and g as 

g, R + (N- 2)g/R = 0, (35) 

and adopt this Robin boundary condition at R = R max . 

We used the finite differencing method with the fourth-order accuracy with respect to 
the grid size. Here, the method of nested hierarchical grids was adopted in the numerical 
calculation. Specifically, we located the outer boundary at -R max = 1024, and put 12 and 16 
layers for < «gb < 10 2 and 10 2 < ogb, respectively (remember that the unit of the length 
is Rs(Mq)). The boundary of the n-th layer is located at R = 1024/2 n_1 , and the grid 
size is 12.8/2 n_1 . Then, the solutions were obtained by using the successive-over-relaxation 
(SOR) method: For each calculation, a test surface is prepared initially and it is made 
slowly converge to the real solution until the difference from the finite difference equations 
normalized by the absolute value of g becomes less than 10~ 12 . 

Figure [T] shows the behavior of g as a function of R for ogb = and 10 fc for D = 5-8, 
where k = 1, 4 for D = 5 and 6 and k = 1, 3 for D = 7 and 8. The curve for «gb = in 
the case D = 5 agrees with Eq. (133]) . As ctGB is increased, the curve becomes steeper around 
the origin. This is the reason why we increased the number of the layers for 100 < q;gb : 
More grid number is required for larger c^b- 

In order to check that the geometry of the generated initial data agrees with that of the 
time-symmetric Cauchy surface in the spherically-symmetric spacetime, Eqs. (130]) and (13T1) . 
we calculate the minimal surface (or the AH) and compare the relations between the two 
nondimensional quantities, oigb I ^(M) and rn/rs(M), where tr is the horizon radius. The 
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FIG. 1: The behavior of g (in the unit of R^(Mq)) as functions of R/ Rs(Mq) for D = 5, 6, 7, and 
8 for the cases a = and 10 k with < k < 3. For D = 5 and 6, the case a = 10 4 is also shown. 

horizon of the analytic solution is given by /(Vff) = 0, which leads to the relation 



r_H_ 

rs 



N-i 



rs 



1. 



(36) 



On the other hand, we find the minimal surface of the numerical data by searching the 
location R = Rh at which 



iV-2" 

is satisfied, and calculate the horizon radius as 



R-qj R + ^ = 



(37) 



r H = R H ^ N ~ 2 \R H ). 



(38) 



Then, the ADM mass M is calculated by Eq. (I26p . and the value of rn/rs{M) is evaluated. 

Figure [2] shows the relation between a G B/r 2 s (M) and r H /r s (M) for D = 5-8. The solid 
curve shows the analytic formula, and the circles are our numerical data points. They agree 
very well, and thus, we can confirm that the generated initial data is the time-symmetric 
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£>=5 D=6 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 1 2 3 4 5 6 7 



a GB /rj(M) a GB /rj(M) 

FIG. 2: The radius rn/rsiM) of the minimal surface as a function of aGB/^gC-^) for D = 5, 6, 
7, and 8. Each straight curve shows the analytic relation, while the circles (o) show the numerical 
result. The numerical data are taken for a^^. / B? s {Mq) = 0, n x lO^" 1 (where n and k are integers 
satisfying 1 < n < 9 and < k < k max ), and 10 fcmax , where k max = 4 for D = 5 and 6 whereas 
&max = 3 for D = 7 and 8. The data agree well with the analytic relation. 

slice in the spherically-symmetric spacetime. The cases for aQ^,/R 2 s {Mo) = 0, n x 10 fe_1 
(where n and k are integers satisfying 1 < n < 9 and < k < k max ), and 10 fcmax are shown 
for all D. Here, /c max is chosen as k max = 4 for D = 5 and 6, and /c max = 3 for D = 7 and 8. 
As seen from this figure, the value of rn/rs(M) is decreased as aGB/Rs(M ) is increased, 
but even at a GB /Rs(M ) = 10 4 , the value of r H /r s (M) is ~ 0.36 and 0.13 for D = 5 and 
6, respectively. Therefore, generating the initial data with a horizon radius r# that is much 
smaller than the Schwarzschild radius requires a very large number of q;gb/-R|(-^o)) an d 
thus, it is a hard task at least in this approach. 

Here, we point out the importance of the location of the outer boundary -R max - Figure [3] 
shows the numerical results with the choices of -R max = 16 and 64. The agreement with 
the analytic relation is fairly good for small «gb, but the deviation becomes significant for 
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FIG. 3: The numerical data taken for -R max = 16 (crosses, x) and 64 (squares, □) in the case of 
D = 5. The deviation from the analytic relation becomes significant for uqb / ^(Mq) > 10 and 
100 for -R max = 16 and 64, respectively, indicating that the larger number of i? max is required for 
accurate numerical calculations as ogb is increased. 

&gb / Rs(M ) > 10 and 100 for -R max = 16 and 64, respectively. This is because the integrand 
of Eq. (|26|) cannot be ignored for R > R max and therefore the ADM mass M is evaluated to 
be a smaller value than the true value. We have to remember that the larger value of -R max 
is required for a larger number of «gb for accurate calculations. 

Now we summarize the lessons obtained in the analysis in this section: (i) The source 
term for the equation g becomes zero at R = because of the cancellation of divergent 
terms, and thus, the regular solution can be constructed; (ii) The ADM mass M has to be 
evaluated after the function g is generated, because it contributes to the mass; (iii) The 
function g becomes steeper around R = as «gb is increased, and therefore, the better 
resolution is required; and (iv) We have to take care of the location of the outer boundary 
-Rmax for large «gb values since the integrand of Eq. (126"|) becomes large at the distant region. 

V. TWO-BLACK-HOLE INITIAL DATA 

Now we turn our attention to the study on the two-black-hole initial data. After gener- 
ating the conformal factor, we study the common AH and assess whether the Penrose-like 
inequalities are held in this system. 
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A. Numerical calculation 



We assume the two black holes to have the same mass and the initial space to be axisym- 
metric with 0(N — 1) symmetry (i.e., there is the z-axis and the directions orthogonal to 
axis have the same structure). To be specific, we introduce the (z, p) coordinates where the 
metric of the flat space is 

ds 2 = dz 2 + dp 2 + p 2 dtt 2 N _ 2 , (39) 

and assume $ = "$?(z, p). The GR solution \l/ is adopted as the Brill-Lindquist solution 
with two equal-mass black holes that is described as 

*o = 1 + \[Rs(M r- 2 (^L. + _L_) . (40) 
Here, the punctures are located at z = ±z Q on the z axis and 

R± ■■= v / M 7 +? ) (41) 

and Rs(Mo) is defined in Eq. ( 129]) . In this setup, the space possesses the two throats and 
three asymptotically flat regions (say, one upper sheet and two lower sheets). The input 
parameters in the numerical calculations are z and «gb i n the unit Rs{M ) = 1. 

In the numerical calculation of g(z,p), we write down the Laplace operator D a D a and 
the source term 5* in terms of {z,p), where the divergent terms of Eqs. (fT4"l) and ffT5]) (i.e., 
s(°) and s^) are canceled out. Since the two black holes have the same mass, there is 
a mirror symmetry with respect to z — 0. For this reason, we choose the computation 
domain as < z < z + A;z max and < p < Ap max where Az max and Ap max are chosen as 
Az max = Ap max = 1024 in the unit Rs(M ) = 1. At the outer boundary, we impose the 
same boundary condition (13^|) as the spherically symmetric case but rewritten in the (z, p) 
coordinates. Similarly to the case of one-black-hole initial data, we used the fourth-order 
finite differencing and the method of nested hierarchical grids. We put 13 layers to the 
computational domain, where the n-th layer has the boundary at z = zq ± Az maiX /2 n ~ 1 and 
p = Ap max /2™ -1 . If the layer crosses z = 0, the region z < is discarded. The grid number 
iVg rid of p coordinate of each layer is varied as 10, 20, and 40 depending on the situation, 
and the grids of z coordinate have the same size. Then, the solution is obtained by the SOR 
method. The relaxation is continued until the deviation from the finite difference equations 
normalized by the absolute value of g becomes less than 10~ 12 . 
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FIG. 4: Relation between numerical error in g(z,p) and inverse of the grid number (l/./Vg r id) 
(proportional to the grid size) for «gb = (x), 1 (□), 10 (•), and 100 (o) for the case D = 5 and 
z = 1. 

The validity of the numerical computation is checked in three manners. First, it is checked 
that the numerical solution g(z,p) in the case of z = agrees with that of one-black-hole 
initial data. Next, as Zq is increased, the solution is expected to asymptote to that of one 
isolated black hole with half mass in the neighborhood of each puncture, and the numerical 
solution g(z, p) is actually confirmed to show this behavior. Finally, we checked whether the 
numerical solution shows the appropriate convergence with respect to the grid size. Figure H] 
shows the typical numerical error in g(z,p) for ogb = 0, 1, 10, and 100 for the case D = 5 
and function of log 10 (l/iVg ric i). Here, the error is evaluated for N gri ^ = 10, 20, 

and 40 by calculating the difference from the solution of the case iV gric i = 80. The slope for 
the curves of «gb = and 100 is approximately four, reflecting the fact that the adopted 
method is the fourth-order accuracy scheme. On the other hand, the slope for the curves 
of a = 1 and 10 is larger than four; the slope for a = 1 is approximately five. This is a 
somewhat strange result, because the five-order accuracy was obtained by using the scheme 
with the fourth-order accuracy. This is probably because the numerical error in g changes 
the value of the right-hand side S of Eq. (1121) . and the change in S further modify the value 
of g. Because of this effect, the cancellation of the numerical error would have happened 
accidentally. Anyway, the numerical data shows at least the fourth-order convergence, and 
this supports the accuracy of our numerical calculation. 

Figure [5] shows the 3D plot of the generated data g(z,p) for z = 1 and q;gb = 0, 1, and 
10 in the case D = 5. Here, the data of the ninth layer (layer of n = 9) are used to draw this 
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Q-GB=0 




FIG. 5: 3D plot of the function g(z, p) in the domain < z < 3 and < p < 2 in the case Zq = 1 
and D = 5. The cases for oqb = (left), 1 (center), and 10 (right) are shown. 

figure. As the value of ckgb is increased, the surface becomes steeper around the puncture. 
B. Common apparent horizon 

The AH is defined as the outermost marginally trapped surface, and it satisfies the 
equation of zero expansion, 9 + = V^k^ = 0, where k^ is a tangent vector of the null geodesic 
congruence from the AH. In GR, the formation of an AH implies the existence of an event 
horizon (EH) outside of it assuming the cosmic censorship and the null energy condition for 
the energy-momentum tensor. On the other hand, in GB gravity, this statement does not 
hold because whether — "H^ obeys the null energy condition is quite uncertain. However, 
also in GB gravity, the formation of an AH at least implies the existence of a region where 
gravity is strong. Furthermore, many theorems, such as the area theorem and the second law 
for a future trapping horizon, has been shown to hold also in GB gravity for a spherically- 
symmetric system with matter of GR branch [49J]. For this reason, the AH may imply the 
black hole formation also in GB in a certain condition. For this reason, it is interesting to 
study the formation of an AH also in GB gravity. 

Because the space is momentarily static in our setup, the equation of an apparent horizon 
is reduced to DiS 1 = 0, where s l is a unit normal to the horizon. Assuming the location of 
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FIG. 6: Coordinate shapes of the AH on a (z,x l ) plane in the unit of Rs(Mq) = 1 in the case of 
D = 5 and ogb = 0, 1, 10, and 100. The cases for several values of z$ are shown with 0.2 intervals 



starting from zero, and the case zq 



(crit) 



is also shown. The location of the punctures are shown 



by dots. As zq is increased, the horizon becomes distorted. 

the AH to be R = h(8), where 9 = arctan(p/z), the horizon equation becomes 

hi 



h ee -(D-2)(h 2 + h 2 e ) 



1 

D-3 V h h 



hi 

+ h e 1 + -tZ 



2(D - 2) 



(D - 3) cot 9 



0. (42) 



h 2 J [(D-3) ^ 

The common AH that encloses the two black holes is found using the so-called shooting 
method by solving this equation under the boundary condition = at 9 = and n/2. 

Figure E] shows examples of coordinate shape of the common AH on the (z, x l ) plane, 
where x l is one of the orthogonal directions to the z axis. Here, the cases D = 5 and 
«gb = 0, 1, 10, and 100 are shown, and the values of zq are 0, 0.2,..., and z^ nt \ where z£ T1 ^ 
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FIG. 7: The regions where the AH can be found (AH formation) and cannot be found (NO AH) 
and contours of M/Mq on a (-zoj^gb) plane in the unit Rs(Mq). The AH cannot be found for 
zo > z^ Tlt \ and numerical data of z^ nt ^ are shown by squares (□). The ADM mass becomes large 
as «gb is increased. 



is the critical value for the AH formation. As the value of zq is increased, the AH becomes 
more distorted. There is the other solution to the AH equation, which corresponds to the 
inner boundary of the trapped region. At z = z^ Tlt \ the two solutions degenerate and the 
solution vanishes for z > z^ nt \ 

Figure [7] shows the region where the AH can be found on the (z , «gb) plane for D = 5- 
8. Here, the unit of the length is adopted as Rs(M ). In general dimensions, the AH is 
not formed for sufficiently large zq. In the cases D — 5 and 6, the value of z { crit) /R S (M ) 
becomes large as c*gb is increased. At first glance, one may think that increasing the coupling 
constant acB helps the AH formation. However, this interpretation is not correct, because 
in this figure, the artificial mass M is used in the length unit Rs{M ). As we have seen 
in Sec. IIII Dl the ADM mass M is changed as «gb is increased following Eq. fl2T)|) . Several 
contours of M/Mq are shown in Fig. [JJ As the value of ogb is increased, the mass M also 
becomes large. 
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FIG. 8: The regions where the AH can be found (AH) and cannot be found (NO AH) on a (zq, ckgb) 
plane in the unit r${M). The numerical data are shown by squares (□), and solid line shows the 
border line of the region of AH formation drawn by interpolation. The dotted line is expected 
border line for large «gb for which numerical data were not taken. The border line is expected to 
intersect zq axis at aQB/f 2 s {M) = 0.5 in the case D = 5, and not to cross the zq axis for D > 6. 

Figure [8] shows the region of AH formation on the (^o ; «gb) plane, but now the 
Schwarzschild radius of the ADM mass, rs(M), is used as the unit of the length. In all 
dimensions D = 5-8, the value z^ nt ^ decreases as the value of ctGB is increased, indicating 
that the coupling constant ckgb makes the AH formation difficult compared to the GR case. 
This is a natural result, because in the spherically symmetric case, the ratio of the horizon 
radius to the Schwarzschild radius, rH(M)/rs{M), decreases as ogb / 'fg(M) is increased. 
The border of the region of the AH formation is shown by a solid curve by interpolating the 
numerical data. The dashed line shows the expected border for large «gb for which numer- 
ical calculation has not been done in this paper. In the case D = 5, it is naturally expected 
that the border line crosses zq axis at aGB/ r I(^) = 0.5, because for this value, the horizon 
radius of the spherically-symmetric black hole becomes zero. For the other dimensions, the 
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border line would not cross zq axis, since r#(M) > for arbitrary value of c*gb (see Fig. |2J). 



C. Penrose inequalities 

The Penrose inequality in GR, 

A A H<^D-2[r s (M)} D -\ (43) 

conjectures that the AH area is bounded from above by the horizon area of a spherically- 
symmetric black hole (Schwarzschild-Tangherlini black hole) with the same mass. The reason 
for this conjecture as follows. If the cosmic censorship holds, an AH is formed in an EH, 
and since the EH is located outside of the AH, its area is expected to be larger than that 
of the AH. Because of the area theorem by Hawking, the area of the EH will increase and 
asymptote to that of a final stationary black hole. Since the horizon area of a rotating black 
hole is smaller than that of a Schwarzschild black hole with the same mass, and the final 
mass is smaller than the ADM mass because of the gravitational radiation, the inequality 
(|43p is expected to hold. On the other hand, if the system with an AH whose area is greater 
than the bound of (1431) . the validity of the cosmic censorship, one of the assumptions of the 
above discussion, may be suspected. 

The Penrose inequality has been attracting^ a lot of attentions, and it was proved for 



momentarily static initial data for D = 4-7 



451 ]. On the other hand, a "counterexample" 



(but not in a strict sense) has been found in Ref. 50j in a system consisting of combined 
portions of the Schwarzschild and Oppenheimer-Snyder spacetimes. This example does not 
contradict the cosmic censorship, but contradicts the assumption that the area of the EH 
is larger than that of the AH. But the Penrose inequality can be modified to match this 
counterexample as follows. The AH in a usual sense can be regarded as a "black hole 
AH" since it is formed in a black hole. On the other hand, we can consider a "white hole 
AH" formed in a white hole whose past-directed outgoing null geodesic congruence has zero 
expansion. In the above counterexample, the black hole AH is located inside of the white 
hole AH, and the area of the white hole AH satisfies the Penrose inequality. Therefore, if we 
adopt the outermost horizon, the Penrose inequality still holds. This Penrose inequality for 
the outermost horizon can be proved in spherically-symmetric case assuming weak energy 



condition 



5l|- 
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In GB gravity, the relation between the Penrose inequality and the cosmic censorship 
becomes unclear because the above reasoning for the Penrose inequality may not hold since 
may violate the energy condition. However, it is still of interest whether the universal 
relation like ( I43p holds or not from a mathematical point of view. The proof for the GR 
case cannot be applied to the GB case at least straightforwardly, because it is unclear if 
—H^u satisfies the energy condition. Therefore, it is interesting to test the inequality using 
the initial data constructed in this paper. 

Other than the original version of the Penrose inequality (j4"5|) . we can consider another 
inequality that may be expected to hold in GB gravity. 1 Namely, since the horizon radius 
r#(M) of a spherically-symmetric black hole is different from rs(M), the Penrose inequality 
may be modified as 

A AH <n D _ 2 [r H (M)} D - 2 } (44) 

where r#(M) is defined by Eq. (1361) . In the following, we test if these two inequalities ( 143]) 
and f l44"|) hold in our system. For this purpose, we define Pi := [ r s{M)] D ~ 2 and 

P 2 : = [ r h{M)} D ~ 2 and evaluate these values for selected values of z and «gb- 

Figure M shows the behavior of Pi as a function of z Q . Here the cases of ogb — 0, 1, 
10, and 100 are shown for D = 5 and 6, and the cases of «gb = 0, 0.2, 2, 20 are shown 
for D = 7 and 8. In all cases, the values of Pi are smaller than unity, suggesting that the 
inequality f )43|) is kept in this system. We find that Pi becomes smaller as q;gb is increased, 
and therefore, the coupling constant tends to help the AH satisfy the inequality ( 143]) if « G b 
is positive. 

Figure [TU] shows the behavior of Pi as a function of z . Here the cases of «gb — 1, 10, 
and 100 are shown for D = 5 and 6, and the cases of «gb = 0.2, 2, 20 are shown for D = 7 
and 8. For all D, the value of P 2 is unity for zq = because the space agrees with the 
time-symmetric slice of the spherically symmetric black hole spacetime. As zq is increased, 
the value of P 2 becomes smaller for all values of a. Therefore, the distortion of the AH 
makes the AH area smaller, and the modified version of the Penrose inequality also holds in 
our system. 

To summarize, the black hole initial data in this paper satisfy both two Penrose in- 
equalities ( l4"3j) and fPf4"]) . and no counterexample has been detected. Therefore, the Penrose 

1 The author thanks Tetsuya Shiromizu for this point. 
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FIG. 9: The behavior of Pi := A AU /n D - 2 [rs(M)} D - 2 as a function of z for a G B = (•), 1 (x), 
10 (o), and 100 (□) in the cases D = 5 and 6 and for a GB = (•), 0.2 (x), 2 (o), and 20 (□) in 
the cases D = 7 and 8. The value of Pi is not greater than unity in all cases, suggesting that the 
Penrose inequality Pi < 1 holds in this system. 



inequalities might hold also in GB gravity under appropriate conditions. 



VI. SUMMARY AND DISCUSSION 



In this paper, we studied the method for generating the initial data for one-black-hole 
and two-black-hole systems in GB gravity. Assuming the initial space to be momentarily 
static and conformally flat, the highly nonlinear equation of the Hamiltonian constraint in 
the N + 1 formalism was successfully solved numerically. Using the generated initial data, 
we studied the common AH that encloses the two black holes, and discussed the Penrose 
inequalities in GB gravity. The result suggests that both two inequalities ( H3|) and ( H4|) hold 
in this system. 

Here, let us discuss whether the proposed conjecture for the condition of the AH forma- 
tion in GR can be generalized to GB gravity. In four- dimensional GR, the hoop conjecture 
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FIG. 10: The behavior of P 2 := A AR /Q D ^ 2 [r H (M)] D - 2 as a function of z for a G B = 1 (x), 10 
(o), and 100 (□) in the cases D = 5 and 6 and for «gb = 0.2 (x), 2 (o), and 20 (□) in the cases 
D = 7 and 8. The value of P 2 is unity for zq = and decreases as zq is increased. Therefore the 
modified version of the Penrose inequality P 2 < 1 also holds in this system. 



52( | is well known as the condition of the horizon formation. The hoop conjecture states that 
a black hole with horizon forms when and only when the hoop length C for a system satisfies 
C < 2irrs{M). This conjecture is loosely formulated in the sense that the definitions of the 
horizon (AH or EH), the hoop, and the mass are not explicitly specified. But this conjecture 
is known to give the approximate condition of the AH formation (under appropriate defini- 
tions of the mass and the hoop which may depend on researchers), although not explicitly 
proved. The point in this conjecture is that the typical one-dimensional length of the system 
is restricted from above if an AH is formed, implying that an arbitrarily long AH does not 
form. But this statement holds only for D — 4, because the black hole can be arbitrarily 
long in higher dimensions as expected from the black string solution, and exp licitly shown 
in [53(. Instead, Ida and Nakao proposed the generalized hoop conjecture |53| as 

C D -3<n D _ 3 [r s (M)} D - 3 , (45) 

where Cd-3 is the typical (D — 3)-dimensional quantity ( "hyperhoop" ) in this system. 
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This hyperhoop conjecture was discussed in several papers [42|, I53M56I] . and the results 
support its effectiveness. Let us consider a two-black-hole system with total mass M, 
and suppose the two black holes to be momentarily at rest with a typical distance L. 
In four dimensions, the typical hoop length is estimated by C ~ 2iirs(M/2) + 2L. In 
a similar manner, in D dimensions, the typical hyperhoop quantity would be Cos — 
n D - 3 [r s (M/2)] D - 3 + tt D - 4 [r s (M/2)] D -*L. Substituting this formula into Eq. flU, one 
has 

This gives an approximate condition for the AH formation at least in a qualitative sense. 

Can the hyperhoop conjecture be further extended to GB gravity? One simple manner 
of generalization would be to change from to in Eq. (145)) as 

C D - 3 <n D ^{r H (M)} D - 3 , (47) 

where r#(M) is defined in Eq. (156% However, this simple generalization does not work 
for a two-black-hole system. Since the value of in this case is given by Cos — 

tt D _ 3 [r H (M/2)} D - 3 + n D ^[r H (M/2)} D ~ 4 L, the condition (jUD gives 

T < (Vd-A [r H (M)] D - 3 ~ MM/2)p- 3 
~ V^-J MM/2)P~4 • I ) 

Let us consider the case D = 5 and suppose ogb and M satisfy the relation o;gb = 
(l/2)r|(M/2). In this case, r H (M/2) = and r H (M) = r s {M)/V2, and Eq. (HE} gives 
L < oo. Therefore, the condition (j4Tj) predicts that the common AH forms even for an arbi- 
trarily long distance between the two black holes. This obviously contradicts our numerical 
result in Fig. [SJ Therefore, the condition for the AH formation cannot be obtained at least 
by a straightforward extension of the hyperhoop conjecture, and a further study is required. 

The numerical work in this paper is the first step toward simulations of black holes in 
GB gravity, and a lot of extensions can be considered. For example, it is necessary to 
extend our method of generating time-symmetric initial data to the method of generating 
time-asymmetric initial data like a boosted black hole. For this purpose, the extension of 
the Bowen-York method 48( should be an interesting possibility. Also, the time evolution 
of black hole initial data could be done using the approximation analogous to close-limit or 
close-slow methods in GR 



43. 



57 



The final goal would be to develop numerically stable 



formulations of numerical GB gravity by extending N + 1 formalism and simulate black hole 
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systems fully numerically to clarify a lot of interesting phenomena such as time evolution of 
an unstable GB black hole, rotating systems, and high- velocity collision of black holes. 
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